Molecular Biology and Evolution

Reproductive isolation can evolve rapidly even among populations with extensive genetic similarity, but the mechanisms underlying this divergence remain an open question in evolutionary biology. Here, we investigate the genetic basis of floral scent divergence, a key trait mediating differential pollinator attraction, in a group of nascent Gongora orchid lineages that occur in sympatry. We combined whole-genome resequencing and floral transcriptomics to characterize the transition from the monoterpene-emitting 'M chemotype' to the pollinator-isolated phenylpropanoid-emitting 'A chemotype' nested within it. While overall genetic differentiation was low, a few highly differentiated regions contained candidate genes in biosynthetic pathways. Of particular interest was a putative o-methyltransferase exhibiting evidence of a recent selective sweep in the A chemotype, with much higher expression in scent-emitting floral tissues of A than in M. Transcriptomic comparison with a more distantly related monoterpene-emitting lineage ('X chemotype') further revealed the coordinated and elevated expression of phenylpropanoid biosynthetic genes in M despite its absence of aromatic floral volatiles, including a putative key enzyme potentially involved in directing metabolic flux towards monolignol biosynthesis. Together, our results support a model in which pre-existing activity in a metabolic pathway may facilitate the evolution of a novel floral scent phenotype, potentially making transitions across the adaptive landscape of scent more accessible than others. We propose that such evolutionary contingency may help explain how divergent scent profiles, and consequently, reproductive isolation, can evolve rapidly via few genetic changes.

Complex traits that span multiple tissues and systems often integrate large numbers of genes across development, physiology, and behavior, making it challenging to identify their essential components. Blood feeding in mosquitoes is one such trait. It is ancestral to the mosquito family, maintained in most species for ~200 million years, and was independently lost in three lineages. These convergent losses offer a natural experiment to discover the genetic, physiological, and neural features required for blood feeding. We assembled high-quality, chromosome-level genomes for seven mosquito species, along with whole-brain tomographic reconstructions. Our study spanned the three known non-blood-feeding lineages (Toxorhynchites rutilus, Topomyia yanbarensis, and Malaya genurostris), blood-feeding relatives, and the variable blood feeder Wyeomyia smithii. Comparing orthologous gene clades, we detected convergent gene loss specific to the three lineages that had lost blood feeding. The losses include the salivary platelet-aggregation inhibitor Aegyptin, blood-activated serine proteases such as Chymotrypsin-1 and 2, and a carboxylesterase expressed in the female fat body and brain glia. The loss of blood feeding also extended to chemosensation. Non-blood feeders lack two odorant-binding protein clades, two ionotropic receptor clades associated with blood-component taste detection, and odorant receptor clades expressed in a discrete, strongly female-biased population of antennal neurons. Female-biased head gene expression was reduced in non-blood feeders. Finally, examination of whole-brain tomographic reconstructions across the species revealed smaller antennal lobes in non-blood-feeding females, consistent with reduced olfactory input. Together, these findings identify a compact set of genes, expression patterns, and brain regions associated with blood feeding, offering an evolutionary entry point for functional dissection of how this complex and dangerous trait is built and dismantled.

The southern house mosquito Culex quinquefasciatus is a major disease vector for multiple arboviruses, including West Nile virus and St. Louis encephalitis virus. Insecticides are widely used to control Cx. quinquefasciatus populations, but resistant mosquitoes have emerged in numerous locations across the globe, including the United States. Much progress has been made in recent years in identifying insecticide resistance (IR) mutations in various pest species. In most species, these studies have focused primarily on SNPs as they are readily detectable. However, several genomic structural variants (SVs) have been experimentally shown to impact IR in Drosophila melanogaster, and structural variation has been observed in natural populations of vector species at a few genes that may be relevant to IR (e.g., cytochrome P450s). Because most population genomic data are derived from short sequencing reads, which fail to detect a large fraction of SVs, the contribution of structural variants to the evolution of IR and other adaptations in vector species is likely to be underestimated in the absence of long-read sequencing data. Here, we present the first genome-wide assessment of structural variants in Cx. quinquefasciatus from wild populations in Puerto Rico and Zambia using PacBio long-read sequencing. We constructed high-quality genome assemblies for each and found a total of 259,046 SVs which consisted predominantly of insertions and deletions. We found that 154,071 of these variants contained or overlapped genes. Of these SVs, 5,225 intersected IR genes and 56 of these were predicted to result in changes to downstream protein function. Several SVs predicted to severely impact protein-coding genes also coincided with signatures of selective sweeps, highlighting their possible role as targets of selection. Many of these variants would have gone undetected using short-read approaches, underscoring the importance of long-read sequencing for comprehensive population genetic analyses of adaptation to pesticides and other strong selective pressures, not only in insects but across the tree of life.

Introgression between highly divergent species is generally selected against, yet some foreign alleles can break through species barriers when they provide a sufficient selective advantage, a process known as adaptive introgression. Here, we characterize the structural variant underlying adaptive introgression from Ciona robusta into C. intestinalis, two species separated by 14.5% net molecular divergence. Previous work identified a narrow peak of C. robusta ancestry on chromosome 5 that reached high frequency in several C. intestinalis populations despite otherwise restricted introgression. Using whole-genome Nanopore sequences from 15 C. intestinalis and three C. robusta individuals, together with short-read datasets including parent-offspring trios, we narrowed down the core introgressed segment from 16 to three genes and characterized copy-number variation affecting the cytochrome P450 gene CYP2U1c. The introgressed haplotype carries a tandem amplification of a 5,173-bp segment containing CYP2U1c, with 23 consecutive copies and conserved repeat boundaries. The nearly 120-kb array could not be fully resolved by long-read sequencing, but coverage consistently supported approximately 23 copies per introgressed haplotype. Despite its extreme size and repetitive structure, the allele showed Mendelian inheritance. CYP2U1c copy number was also polymorphic in C. robusta. Long reads fully resolved two structurally distinct haplotypes in one individual, carrying either a duplication encompassing adjacent CYP2U1b and CYP2U1c genes or five CYP2U1c copies. Their distinct breakpoints indicate independent origins. Our results reveal independently arising CYP2U1c amplifications in the donor species and show that an extreme high-copy-number allele was subsequently favoured following introgression into C. intestinalis.

Universal Single-Copy Orthologs (USCOs) are taxonomically-defined gene sets originally proposed to evaluate completeness in genome assemblies. Due to their single- or low-copy nature, they have also been proposed as ideal markers for phylogenomics. Here, we test this assumption using two vertebrate datasets, mammals and ray-finned fish, and their corresponding USCO reference marker sets, actinopterygii_odb10 and mammalia_odb10. Our results show that although most USCO markers are single-copy, 85% and 91% respectively, empirical analyses using a large number of taxa are able to identify a non-negligible proportion of paralogs in the actinopterygii_odb10 and mammalia_odb10 BUSCO gene sets. Gene phylogenies show that several of these paralogs derive from ancient whole-genome duplication events, particularly the teleost-specific whole-genome duplication, and also more recent duplications such as those in salmonids and cyprinids. Paralogs do not reflect the organismal phylogeny but rather the history of gene duplications, so retaining paralogs in phylogenomic datasets can lead to inaccurate phylogenomic hypotheses. To test the effect of ignoring paralogy in BUSCO gene sets, we analyze the mammal and ray-finned fish datasets using concatenated maximum likelihood and summary coalescent methods by randomly choosing one paralog per locus if these are present. The resulting phylogenomic trees are overall consistent with the reference phylogenies, although small discrepancies were found. Phylogeny-guided paralogy removal also changes the sequences retained for many loci, but whether these changes alter the inferred species trees remains to be tested.

Signals of positive selection around the SLC30A9 gene have been reported in human populations outside Africa. Selection likely acted on a highly differentiated single-nucleotide polymorphism, rs1047626, leading to a non-synonymous substitution in the encoded zinc transporter. Because of the striking similarity between the putatively selected SLC30A9 haplotype observed in several current human populations and the Denisovan individual, previous work has proposed adaptive introgression. Yet alternative explanations, including ancient human variation, and the precise archaic source -Neanderthal or Denisovan- remained unresolved. Considering the potentially complex evolution of SLC30A9, we applied Approximate Bayesian Computation (ABC) algorithms coupled to machine learning to investigate the most plausible evolutionary origin of this substitution. After modelling different evolutionary scenarios with forward-in-time simulations, our results highlight that the most probable scenario is a Denisovan origin of the rs1047626 polymorphism. However, the allele likely introgressed into Neanderthals first and was then passed into non-African modern humans. Moreover, the derived allele frequency for rs1047626 across several African populations is consistent with back-to-Africa migrations. Finally, our ABC analyses indicate strong positive selection in East Asian populations and other out-of-Africa populations, whereas in Neanderthal populations, the selection coefficient was probably neutral or slightly deleterious.

Understanding how genomic architecture and demographic history interact to shape patterns of genomic divergence is central to explaining adaptive radiation and speciation. Here, we provide the first full sequence patterns of diversity and differentiation for Hydrophis sea snakes, using population sampling of three ecomorphologically distinct, coexisting species that diverged near-simultaneously within the last few million years (H. major, H. elegans and H. curtus). We found recent population expansions in all three species, with little evidence for pervasive interspecific gene flow. Genomic landscapes of differentiation were correlated across pairwise comparisons, with relative divergence associated with local reductions in nucleotide diversity but mostly decoupled from absolute divergence. Peaks of differentiation were concentrated in putative pericentromeric regions, consistent with linked selection in shared low-recombination regions. In contrast, a large region of low relative differentiation on chromosome 2 showed elevated diversity across all species and overlapped a prominent cluster of chemosensory vomeronasal type-2 receptor (V2R) genes. Overall, our results suggest that shared demographic history and linked selection shape a common genomic landscape in Hydrophis, while chemosensory receptor gene evolution may have contributed to ecological diversification during the land-to-sea transition of these unique marine reptiles.

Whole-genome duplications (WGDs) have long been proposed as evolutionary facilitators of biological complexity and diversification. Ancient WGDs occurred deep in the evolutionary history of both plants and vertebrates, but appear to be rare in invertebrates. However, it is unclear whether this reflects genuine scarcity or inadequate detection methodologies, and current methods may be insufficient to detect very old WGDs where extensive gene loss has followed. Here, we used a combination of paralogous and orthologous synteny-based methods incorporating bilaterian ancestral linkage groups (ALGs) to search for genomic signatures of WGD in members of the invertebrate phylum Bryozoa. This revealed strong evidence for an ancient WGD in the freshwater bryozoan Cristatella mucedo, Class Phylactolaemata, with duplicated paralogous regions across every chromosome. Only ~10% of duplicated genes are retained, including a Hox cluster duplication. All eight genomes from the bryozoan class Gymnolaemata also show signatures of ancient WGD in the form of 2:1 ratios in orthologous synteny comparisons, which remain detectable despite only ~5% duplicate retention. Gene trees support a scenario in which all extant bryozoans share an ancient WGD and much of the genome rediploidised independently in the two clades. Duplicates retained after bryozoan WGD are enriched in cilia genes and are preferentially expressed in the lophophore, the tentacular feeding organ covered in highly specialised cilia. We argue that cryptic WGD in bryozoan evolution facilitated adaptation to a sessile, filter-feeding lifestyle, and propose that incorporation of ALGs adds to the power of synteny-based WGD-detection methods.

Tandemly repeated sequences are noncoding, highly variable, yet integral components of eukaryotic genomes. Empirical research has so far largely focused on a few well-characterized, functional tandem repeats (TRs) or documented the broad diversity and rapid evolution of these repeats without explicitly testing the forces driving their evolution. By characterizing TR abundance and TR sequence composition within and between chromosomes and species in the stick insect genus Timema, we uncover that TR loads are enriched in genomic regions of high recombination but tend to be reduced in species with higher effective population sizes. These results suggest that there is a mutational effect of recombination that generates increased TR loads and that the bulk of TRs experience weak selection, leading to their rapid turnover over short evolutionary timescales. Finally, comparisons of centromere and non-centromere TRs reveal that the rapid divergence of centromeres is best explained by the combination of the fast-evolving nature of TRs per se and centromere-specific evolutionary forces. Overall, this study provides the first in-depth investigation of genome-wide TR evolution, while integrating the effects of recombination, genetic drift, and centromeric position along the genome.

Phenotypic plasticity is a common adaptation to heterogeneous environments. While it is well documented that trait plasticity can rapidly evolve, little is known about the developmental mechanisms that underlie variation in plasticity. Here, we compare diverged populations of the common buckeye butterfly (Junonia coenia): one population with canalized wing patterns that are insensitive to environmental conditions and another population that has distinctive seasonal color patterns. In common garden trials we were surprised to find dramatically elevated levels of transcriptional and chromatin accessibility plasticity during wing development of the canalized population compared to the seasonally plastic population. Further, using machine learning analysis of accessible chromatin sequences, we observed a broad heterochronic shift in motifs enriched in plastic putative regulatory elements to earlier stages of development in the canalized population, which we speculate may cause phenotype-related gene expression to move outside a wing color critical period. Together, our data lead us to propose a model where cryptic genetic processes are masked by a shift in chromatin remodeling, which in turn could produce a capacitor-like effect to maintain developmental genetic variation in the canalized population.

Populations adapting to heterogeneous environments often harbor many alleles that covary with environmental gradients, yet how so many adaptive variants can be maintained and show signatures of local adaptation remains unclear. Epistatic non-additive interactions among alleles may concentrate fitness effects at specific allelic combinations to shape complex traits and facilitate adaptation. We evaluate compositional epistasis among three intronic SNPs in the couch potato (cpo) gene, each independently clinal along the North American latitudinal gradient, using recombinant outbred populations fixed for distinct three-locus combinations. Flies carrying all three temperate alleles (cpoTTA) or all three tropical alleles (cpoCGT) show non-additive increases in traits associated with fitness in their respective environments: dormancy, starvation resistance, and lipid storage in temperate flies and fecundity in tropical flies. Intermediate genotype combinations generally show weaker and negative deviations from additive predictions. Temperate genotypes show reduced cpo expression and independently create predicted transcription factor binding sites, suggesting a regulatory basis for the observed phenotypic effects. Consistent with selection shaping haplotype frequencies against a backdrop of ongoing recombination, the temperate and tropical genotypes are enriched above random expectation across the latitudinal cline, with directionally consistent but non-significant seasonal patterns. Our findings demonstrate that non-additive interactions among non-coding variants within a single gene can produce large, pleiotropic fitness effects and shape genotype distributions in natural populations, providing a mechanism by which a small number of regulatory variants can underlie complex climatic adaptation.

Hox proteins are evolutionary conserved transcription factors guiding morphogenesis in all bilaterian animals. Conservation mainly concerns the homeodomain (HD), a helix-turn-helix DNA binding motif and flanking sequences close by. Besides, sequences are poorly conserved and disordered. Here we investigate how poorly conserved Intrinsically Disordered Regions (IDRs) influence Hox protein activity, focusing on a closely related Drosophila Hox pair and using a wide range of in vivo function. Studying Drosophila Ultrabithorax (Ubx) and AbdominalA (AbdA) proteins, we find that conserved protein motifs close to the HD are differently used in Ubx and AbdA. Coarse-grained molecular dynamics simulations identify a large set of intra-molecular interactions, notably involving poorly conserved disordered regions upstream and downstream the HD. Swapping these regions in Ubx and AbdA alters protein functions and switches motif dependency in some biological situations. This suggests that IDRs, predicted to mediate large and dynamic sets of intramolecular interactions, influence protein motif usage in a context dependent manner.

Phylogenetic trees inferred from different genomic regions, or from organellar and nuclear genomes, frequently disagree owing to incomplete lineage sorting, hybridization, introgression or horizontal gene transfer. Such discordance is informative, yet the standard tool for comparing trees visually, the tanglegram, is restricted to two trees and is known to misrepresent their similarity. Here we introduce phylogenetic parallelograms, drawings in which any number of rooted trees are embedded into a common scaffold, a rooted phylogenetic network that displays all of them, and are drawn in parallel so that shared branches align while conflicting branches visibly diverge. On synthetic tree pairs, the complexity of a parallelogram tracks the topological distance between the trees, whereas the complexity of an optimized tanglegram does not. We apply the approach to genome-wide introgression in the Anopheles gambiae species complex and to organellar discordance in cats, grasses and Fagaceae. PhyloParallelograms is an open-source interactive application.

GC-biased gene conversion (gBGC) is a widespread recombination-associated process that favours the fixation of GC alleles independently of their fitness. In coding regions, it can oppose purifying selection and favour the accumulation of deleterious AT to GC mutations. In mammals where recombination hotspots are predominantly determined by PRDM9, the deleterious effect of gBGC is expected to be episodic due to the rapid turnover of hotspots. The deleterious impact of gBGC episodes on proteins has been well documented, particularly in rodents and primates. However, the evolutionary response of the affected proteins after the episodes end has been scarcely investigated. Theory predicts that a decline in gBGC intensity should trigger a compensatory response that enables previous damage to be repaired. Recently, we reported evidence for pervasive compensatory evolution following gBGC episodes in Hydromyini rodents. Here, we aimed to investigate whether this post- gBGC compensation is widespread among mammals in which gBGC is episodic. Using a phylogenetic approach, we analysed substitution patterns in coding exons from five mammalian taxa spanning a range of effective population size (Ne): Murini, Rattini, Cetacea, Cercopithecidae and Hominoidae. We detected thousands of gBGC episodes, characterised by an increased rate of non-synonymous AT to GC substitutions. We found evidence for compensatory evolution in non-primate taxa, and a lack of signal in Primates. Although the influence of Ne remains unclear, our study highlights both the evolutionary significance of post-gBGC compensation and the methodological challenges in detecting it.

The 21st century has seen repeated challenges to the premise that de novo mutations occur independently from their fitness effects. So far, these attempts have been dismissed as technical artifacts, restoring confidence in the null hypothesis. Yet the development of high-resolution genomic resources enables investigating the relationship between mutation rate variation and selective constraint with unprecedented detail. Here, we use state-of-the art genomic maps to test whether predicted mutation rates differ within and outside putatively constrained sites in humans. Contrary to the hypothesis that is typically put forward, but in agreement with another recent study, we find that protein-coding genes exhibit higher rates of mutation than the neutral genomic background. Conversely, mutation rates are lower in non-coding elements identified through sequence conservation. This can be attributed either to biases in the estimation of mutation rates in sites under selection, biases in scoring constraint in sites with low rates of mutation, or a genuine signal that would merit further investigation. Regardless of the underlying cause, our findings have broad implications for evolutionary genetics. Several studies rely on accurate mutation maps and identification of selection targets across the genome, which inherently confound each other. Caution is therefore warranted when interpreting results from currently available non-coding scores.

Pinpointing the genetic elements underlying a phenotype within a supergene can be difficult because hundreds of linked genes can show similar associations with phenotype. Rare recombination between supergene haplotypes can help identify those genetic elements by changing haplotype-phenotype combinations. Here we investigate the genetic underpinnings of colony queen number and body size in the ant Formica lemani. Throughout the Formica genus, single-queen colonies typically carry M haplotypes, while multi-queen colonies also carry P haplotypes. In F. lemani, we found two P haplotypes homologous to those described in the closely related F. cinerea, but unlike in F. cinerea, where P1 is restricted to multi-queen colonies, F. lemani P1 is common in both single- and multi-queen colonies. This deviation suggests that P1 may carry novel SNP combinations allowing single queen colony compatibility. By comparing haplotypes within F. lemani and with Formica cinerea, we searched for candidate monogyny variants resulting from gene flux from the M onto the P1 haplotype, or with P1-specific mutations. We found 82 genes associated with candidate variants, of which two are among the fourteen previously suggested candidates in other species, suggesting that several different genes may underlie queen number across the genus. We also identified a supergene associated with smaller queens and males, which is shared with F. cinerea and in interchromosomal linkage disequilibrium with the queen-number supergene. Together, these results provide a rare example of how an ancient supergene can change in function, and a step toward identifying the genes within a supergene responsible for its functions.

Insertion sequences (IS) are the most common type of transposable element in prokaryotes and shape the structure of genomes through transposition and by providing a substrate for recombination. Despite the mutational impact of IS, the evolutionary dynamics of most elements in host species remain unknown. Here we study the dynamics of IS6110 in 10,000 strains of the Mycobacterium tuberculosis complex (MTBC). We developed a tool that allows the detection and comparison of IS insertions from short reads without using a reference genome. Using ancestral state reconstruction (ASR) on presence-absence patterns of IS6110, we describe the distribution of copy numbers (CNs) in the MTBC, infer birth rates of the element, and identify genomic regions with large numbers of parallel IS6110 insertions. Copy numbers in the MTBC range from 1 in some clades to more than 30 in strains of La3 (M. orygis). IS6110 birth rates scale approximately linearly with copy number and are elevated on terminal branches, consistent with the delayed action of purifying selection. A key characteristic of IS6110 is its occurrence in hotspots: the 5% most frequently targeted regions account for half of all independent insertion events. The motif 5'-TCTCAAAW-3' is enriched around target sites and in hotspots, suggesting that the accumulation of insertions in these regions results through a combination of non-random insertion and purifying selection in other regions. To conclude the study, we propose a niche constraints model according to which the distribution of IS6110 in the MTBC is governed by the rarity of regions that have both suitable DNA properties and little functional value for the host.

Emergence and evolution of functional RNA and protein structures are the central problems for understanding the origin of life. Although it is well known that catalytically active RNA elements, ribozymes, can catalyze many reactions, including peptide bond formation, the specifics of the transition from the hypothetical, primordial RNA world to protein-based life centered at the translation system remain enigmatic. We developed AMES, Atomistic Molecular Evolution Simulator, and employed it to perform computer simulations of the evolution of short RNA molecules and RNA-peptide complexes. Comparison of the evolutionary trajectories and the structures of RNA molecules and RNA-peptide complexes emerging in these simulations shows that short random peptides accelerate RNA evolution, stabilize RNA folds, and boost the structural diversity of evolving RNA molecules, potentially enabling a broader range of activities. We hypothesize that the RNA world was, actually, an RNA-peptide world, in which, from the earliest stages of evolution, evolving RNA molecules interacted with short random peptides synthesized in a non-templated manner. These interactions could drive the evolution of diverse RNA structures and activities, and eventually, of the translation machinery.

Crozier's Paradox states that genetic kin recognition will select for its own demise, and theoretical analyses over the past few decades have supported this conclusion. However, empirical evidence is far more limited. Here we examine support for Crozier's Paradox in two kin recognition genes in the social amoeba Dictyostelium discoideum, tgrB1 and tgrC1, which enable cells to recognize and reject unrelated cells during cooperative multicellular development. We find that co-isolated amoebae, taken from just a few grams of soil, harbor highly diverged recognition alleles, indicating plentiful sequence variation for social decision-making on a micro-scale. tgrB1 and tgrC1 belong to a large gene family, whose members vary in their levels of polymorphism, strength of selection, and retention of polymorphism. We do not see strong evidence for the paradox, in the form of positive frequency-dependent selection, or its outcome, in the form of depleted genetic variation. Consistent with previous findings, we detect recombination among the duplicated tgrB and tgrC genes, which suggests a mechanism to generate novel alleles and replenish genetic variation. We suggest that the evolutionary dynamics of large gene families may be essential for rapidly generating novel variation that can prevent Crozier's Paradox from taking effect.

Mutation-selection balance is the most plausible general mechanism for the maintenance of additive genetic variance (VA) for quantitative traits in natural populations. It is well established that the genetic variants affecting traits (QTLs) are widely distributed across species' genomes. At the same time, unconditionally deleterious mutations are also expected to be common, with those with the largest fitness effects tending to be strongly recessive. Many genomes possess regions of restricted recombination that could result in tight linkage between largely recessive deleterious alleles, conditions which can lead to elevated heterozygosity via associative overdominance (AOD). Here we show that weakly selected QTLs linked to genomic regions exhibiting AOD can elevate VA above what is expected in the absence of unconditionally deleterious mutations. Selection on traits when VA is maintained by AOD, however, may expose recessive deleterious alleles as homozygotes, reducing fitness. Thus, while AOD may contribute to the maintenance of VA, variation maintained in this way may not be readily accessible by natural or artificial selection.